      SUBROUTINE BFS(L,X,EPS,F,FP,FPP)
c
c23456789012345678901234567890123456789012345678901234567890123456789012
c
C    BFS CALCULATES UNNORMALIZED SPHERICAL BESSEL FUNCTIONS OF THE FIRST
C    KIND( J FUNCTIONS). RECURSION DOWNWARD IN L IS USED. OUTPUT IS THE
C    FUNCTION, F, FIRST DERIVATIVE, FP AND SECOND DERIVATIVE, FPP.
c
c    calls no other routines
c
      IMPLICIT REAL*8(A-H,O-Z)
c
      FL=L
      EM=14.D0
      IF(X.GT.FL) EM=EM+X-FL
      M=EM
      FL1=FL+1.D0
      FL2=FL+FL1+2*M
      F=0.D0
      F1=EPS
      DO 1 I=1,M
      F3=F
      F=F1
      FL2=FL2-2.D0
    1 F1=FL2*F/X-F3
      FP=FL*F/X-F3
      FPP=-2.D0*FP/X+(FL*FL1/(X*X)-1.D0)*F
c
      RETURN
      END
